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SUMMARY 

The methods of the calculus of variations are used to maximize payload capability for 
multistage launch vehicles. The method of solution uses the Lagrange multipliers to de- 
termine the optimum thrust direction profile, as well as to construct partial derivatives 
of payload with respect to the stage propellant loadings and a booster steering parameter. 
These derivatives are used to terminate the stages and/or as terminal equations to be 
satisfied. Maximum payload can thus be achieved with a single solution, rather than 
with a family of parametric results. Constant thrust and specific -impulse operation is 
assumed for each upper stage (booster thrust and specific impulse vary with atmospheric 
pressure), and structure weight can be either fixed or a linear function of the stage pro- 
pellant loading. Two-dimensional flight in a central, inverse -square gravitational field 
is assumed. 

Numerical results are presented for two- and three-stage launch vehicles flown to 
circular orbit and Earth escape, respectively. Parametric results are presented and 
compared with the overall optimum solution obtained by use of the variational technique. 


INTRODUCTION 

A problem that frequently arises in trajectory optimization studies is that of deter- 
mining the maximum payload capability of a multistage launch vehicle flown to a pre- 
scribed set of burnout conditions. If all vehicle parameters are specified, the problem 
reduces to that of finding the optimum steering profile. In many cases, however, not all 
these parameters are specified, and those left unspecified can be varied to maximize 
payload. 

A typical situation that occurs in the design of future launch vehicles is one in which 
the propulsion system (thrust and propellant flow rate) is specified, but some or all the 
stage propellant loadings are left unspecified. The unspecified propellant loadings gen- 
erally can be varied to achieve maximum payload capability for the vehicle. 



An additional optimizing parameter frequently is available in the booster steering 
program. Since the booster stage operates in the atmosphere, the booster thrust direc- 
tion profile is shaped to minimize aerodynamic heating and loads and is not available for 
complete optimization. A single degree of freedom remains, however, corresponding to 
the magnitude of a short pitchover phase following the initial vertical rise. This degree 
of freedom, sometimes called the booster kick angle, determines the amount of trajectory 
lofting during boost phase. Since the upper stages operate essentially under vacuum con- 
ditions, the steering program for these stages is available for complete optimization. 

Many authors (e.g. , refs. 1 to 7) have treated the problem of optimizing the stage 
propellant loadings of multistage vehicles. None of these authors, however, has at- 
tempted to optimize the steering program for these vehicles. Others (e.g. , refs. 8 
to 10) have used the calculus of variations to optimize the steering program for various 
rocket vehicles. In particular, reference 11 treats the problem of optimizing the steer- 
ing program of a multistage launch vehicle. Reference 11, however, does not consider 
the problem of optimizing the stage propellant loadings or booster kick angle. 

Recently, Mason, Dickerson, and Smith (ref. 12) have considered the problem of 
simultaneously optimizing the steering program and the stage propellant loadings of a 
multistage launch vehicle. These authors followed the approach of Denbow (ref. 13) and 
Hunt and Andrus (ref. 14) in formulating the variational problem. 

The present report was written concurrently with reference 12 and presents a method 
which allows the propellant loadings, booster kick angle, and upper -stage steering pro- 
gram to be simultaneously optimized. The variational approach is somewhat different 
from that used in reference 12. By following the method of reference 15, the maximizing 
functional is written as the sum of the final payload and a constraint integral for each of 
the upper stages. The resulting boundary equations supply partial derivatives of payload 
with respect to the unspecified parameters. These derivatives are then used, along with 
the required burnout conditions, as terminal equations to be satisfied. The analysis does 
not require that all the stage propellant loadings (or booster kick angle) be optimized. 
Equations are developed for optimizing payload with respect to any combination of un- 
specified parameters. 

The variational equations for optimizing vehicle parameters have been incorporated 
into a digital computer program used previously at Lewis for parametric launch -vehicle 
studies. This computer program is not discussed in detail; however, numerical results 
are presented for two- and three-stage launch vehicles flown to circular orbit and Earth 
escape, respectively, to demonstrate the feasibility of the variational approach. Para- 
metric results are presented showing the variation of payload with propellant loadings 
and booster kick angle. The resulting payload envelopes are then compared with the over- 
all optimum points generated directly by means of the variational technique in order to 
verify the equations. Along with the results, the procedures used to obtain numerical 
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results are briefly discussed. 


ANALYSIS 

The problem to be solved is to determine the maximum payload capability of an 
N-stage launch vehicle flown to a specified set of burnout conditions. The analysis admits 
atmospheric effects during booster phase but assumes vacuum operation for all other 
stages. Because of these atmospheric effects, the booster steering program is assumed 
completely specified (e.g. , zero angle of attack), except for the booster kick angle. The 
upper-stage steering program, however, is unconstrained and is determined to maximize 
payload. The calculus of variations is used for this purpose. 

Each of the upper stages is assumed to operate at a fixed (constant) value of thrust 
and propellant flow rate. These values may be zero, so that coast phases are admitted. 
The structure mass for each stage is assumed to be a function of the stage propellant 
loading, defined by 


m 


s 


= m H + km p 


where m g is the total structure mass, m^ is the fixed mass, m p is the stage propel- 
lant mass, and k is the propellant sensitive mass fraction. (All symbols are defined in 
appendix A . ) 

In addition to the variational trajectory, provision is also made for adding an addi- 
tional velocity increment Avj after the desired orbit conditions are achieved. This ve- 
locity increment is achieved by use of the final stage for propulsion. The amount of pro- 
pellant required for this maneuver is calculated by use of the standard impulsive velocity 
equations . 


Variational Problem 

Since the booster steering program is not subject to complete optimization, the 
booster stage is not treated in the following Euler -Lagrange equations. The booster de- 
grees of freedom (propellant loading and kick angle) are included by allowing variations 
in the position and velocity at second-stage ignition. The associated equations, along 
with the equations for optimizing upper -stage propellant loadings, are treated in the 
boundary equations resulting from the variational analysis. 

The variational problem to be solved is that of finding the upper -stage thrust pro- 
gram which maximizes the payload capability of an N-stage launch vehicle for given 
boundary conditions. This problem can be formulated as a generalized Bolza problem. 
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By following the treatment in reference 15, 
the functional to be minimized can be 
written as 



a) 


where g is a function of initial and final 
conditions to be minimized and F. con- 
sists of a set of constraints to be applied to 
each of the upper stages, added together 
with the aid of Lagrange multipliers. For 
this problem, the payload is to be maxi- 
mized, so that 


g = -m pL 


( 2 ) 


The constraint equations are 


f 


li 



co 2 r 


T. 

— sin i/j = 0 
m 


(3a) 


T. 

f = rcb + 2uc o — cos ijj = 0 

Z1 m 


applicable on the interval 


f 3i = r - u = 0 


f 4i = <P ~ w = 0 


i 5i =m + ^ = 0 


t. 1 < t < t. 
l-l — — i 


for i = 2, . . . , N 


(3b) 

(3c) 

(3d) 

(3e) 
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Equations (3a) to (3d) are the two-body equations of motion written in two-dimensional 
polar coordinates with an inverse square force field acting. The thrust direction ip and 
the state variables are defined in figure 1. Equation (3e) defines the propellant flow 
rate. 

Equations (3) are combined to give 


F i = E Vii 

3=1 


i = 2. 


N 


(4) 


where Aj ■ are undetermined Lagrange multipliers, which are functions of time since 
the constraint equations must be satisfied at all points of the trajectory. 


Euler-Lagrange Equations 


As shown in reference 16, a necessary condition for g to be minimized is that the 
Euler-Lagrange equations be satisfied. The Euler-Lagrange equations are 


d_ 

dt 



9F i 


9Xj 


i = 2, 


N; j = l, 


(5) 


where 


X j 


are the problem variables 


X x (t) = u 

X 2 (t) = CO 

x 3 (t) = r 
x 4 (t) = <p 
x 5 (t) = m 

Xg(t) = ip 


The Euler-Lagrange equations for the present problem can be written explicitly by use 
of equations (3) and (4): 
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A ^ “ 2 Cl) A 2 " A.0 

(7a) 

* n ^4 

A n — 2Cl)A-J + — Art “ — 

(7b) 

r r 

x 3 - - <Ah + wa 2 

(7c) 

\r / 

X 4 = 0 

(7d) 

. rp 

X 5 ” (*i sin ^ + ^2 cos ^ 

2 

(7e) 

m 

T 

(A^ cos ip - ^2 sin ^) — = 0 

(7f) 


m 


Equations (7) (and subsequent equations) apply separately to each of the upper stages. 
The subscript i has been omitted for simplicity. 

Equation (7f) determines the thrust direction (for T + 0): 


tan 



(8a) 


sin i p = ± 


cos 4>= ± 


A 1 

V X 1 + x 2 
A 2 

+ x 2 


(8b) 


(8c) 


The uncertainty in sign in equations (8b) and (8c) corresponds to an equivalent uncertainty 
of 180° in the thrust direction. The choice of sign will be determined later in the anal- 
ysis. 

Equations (7a) to (7d) must be integrated, along with the equations of motion 
(eqs. (3a) to (3d)) to determine the thrust direction and the optimum trajectory. It is 
shown later that equation (7e) need not be integrated. Equation (7d) is easily integrated 
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to give 


A^ = Constant (9) 

First Integral 

Since the function F does not contain the independent variable (time) explicitly, a 
first integral to the Euler -Lagrange equations exists (for each stage), which can be 
stated as (ref. 16) 


6 


F - 



X i 


= C 


(10) 


Equation (10) can be written explicitly by use of equations (3) and (4): 

C-/JL-J 

V 2 

Equation (11) is used later in the analysis to eliminate both C and Ag from the boundary 
equations. 


r\ Aj - 2 ucoA 2 + uAg + coA^ - /3Ag + — (Aj sin ip + Ag cos ip) = 0 (11) 

/ 


Weierstrass Condition 

The uncertainty in sign in equations (8b) and (8c) can be resolved by applying the 
necessary condition of Weierstrass (ref. 16). Following the development in refer- 
ence 17, this condition can be stated as E > 0 for a minimum, where 

6 

E = F( X * ,i* ) - Flxj.ij) (ij - ij) || (12) 

j=l 3 

The Xj correspond to the minimizing values, which differ from the x? by a finite but 
admissible amount. For the present problem, only tp is subject to such a variation. 
Equation (12) can be evaluated to give 
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Equations (8b) and (8c) are now used along with 


sin jp* = t 




V3 


2 , 2 
+ Xr> 




A o 

cos Ip* = T ■ ~ 


VS 


+ Xr 


( 14 ) 


to give equation (13) in the form 


2T 


m 


+ Ao > o 


2 - 


(15) 


Equation (15) implies that the plus sign must be chosen in equations (8b) and (8c). 


sin ip 


M + 4 


(8b 1) 


cos ip- 


& 


+ A‘ 


(8c 1) 


Weierstrass-Erdmann Corner Condition 

The boundary conditions on the Lagrange multipliers at staging points can be derived 
with the aid of the Weierstrass-Erdmann corner condition (ref. 16). This condition 
states that 0F/3Xj (j = 1, . . . , 6) must be continuous at such corners. For the present 
problem, 



3F 

9u 


= H 


3F 


= rX 


ad) 


2 


3F 

af 


- X 3 


3F 

d<p 


> 


X 4 


3F 

am 


X 5 



dip 


J 


( 16 ) 


All the multipliers are thus continuous across staging, hence continuous throughout the 
trajectory. Equations (16) also imply the continuity of the thrust direction and rate ip 
and ip since these variables are functions of the state conditions and of the Lagrange 
multipliers (appendix B); however, the constant of integration C is not, in general, 
continuous. From equations (11) and (16) the discontinuity in C is given by 


AC = - 





A/3 


(17) 


Transversal ity Equation 

The relation between changes in boundary conditions and changes in J is expressed 
by the general transversality equation (ref. 16). For this problem, the transversality 
equation can be written 
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This equation can be written explicitly by using the definition of F (eqs. (3) and (4)) and 
the first integral (eq. (10)). 



The subscript i has been used with C and dm since these variables may be discon- 
tinuous at staging points. 


Boundary Equations 

If some of the problem variables (state conditions or control variables) are not 
specified, values should be chosen which minimize J (or, equivalently, maximize pay- 
load). 

According to reference 16, minimizing J is accomplished by setting dJ (eq. (19)) 
equal to zero. Equation (19) has the form 

m 

dJ = y Gi dx. = 0 (20) 

j=l 

If the m problem variables Xj are all independent, dJ will vanish if, and only if, each 
term on the right side of equation (20) is independently set equal to zero. For specified 
variables Xj, the allowable variation dxj is zero; for unspecified x^ the coefficient 
Gj must be set equal to zero. Equation (20) can thus be interpreted as a total differential 
of J, and 


G. = °± (21) 

3 9x j 

Equation (19) is not suitable for this interpretation, since the variables are not all inde- 
pendent. In the following section, the dependent variables in equation (19) are eliminated 
by expressing the dependence explicitly. 

Consider first the terms in equation (19) involving the variations of the state vari- 
ables. 
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N t N-l 

f\ 1 du + rx 2 dw + X 3 dr + X 4 d cp) 1 = - ^ Rx-j^ du) + - (Xj du)“ + (rX 2 dw) + - (rX 2 do>) 

i=2 V h-l i=2 L 


(X 3 dr ) + - (x 3 dr)" + (x 4 dcpf - (x 4 dcp ) ^ 


+ (x^ du + rX 2 dw + X 3 dr + X 4 dcp) _ - (X^ du + rX 2 do> + Xg dr + X 4 dcp) + 

t=t N t=t l 


( 22 ) 


where the superscripts - and + refer to conditions before and after staging, respec- 
tively. The state variables are continuous throughout the trajectory, so that 

( dx l) t=t J- ’ Mt.t- 1 = 2 - 1; j = 1, . . . , 4 

Since the Lagrange multipliers are also continuous (eqs. (16)), the terms before and 
after staging in equation (22) are equal and cancel. 

At t = tj, the variations in state conditions are due to the allowable variations in 
booster burning time and kick angle. 


Sx. 3x. 

(dx.) = — 3. dT 1 + — - da j = 1, . . . , 4 (23) 

v ] 't=q dr 1 1 9a 

where a is the booster kick angle and Tj is the booster burning time. (In general, i\ 
is the burning time of stage i, that is, Tj = ^ - t.^. ) The partial derivatives in equa- 
tion (23) are evaluated by using a numerical approximation method, which is discussed 
later . 

The state variations at t = t^, as expressed in equation (22), may or may not be in- 
dependent, depending on the specification of the desired burnout conditions. Some typical 
burnout requirements are presented in appendix C. For generality, the variables Xj 
are expressed in terms of a set of generalized (independent) state variables, r}^, 
k = 1, . . . , 4, so that 


4 



k=l 


0 Xj 


dr k 


d7 k 


j = 1, 


4 


(24) 
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Equations (23) and (24) are combined with equation (22) to give 



(25) 

The superscripts + and - are omitted in equation (25) and subsequent equations. The 
expressions to be evaluated at t^ and t^ should be evaluated at tj and t^, respec- 
tively . 

Inasmuch as a suitable form for the state variation terms in equation (19) has been 
obtained, the mass (and payload) variation terms are now considered. Since the pro- 
pellant flow rate for each stage is constant, the variables rm are not independent and 
can be expressed as functions of the burning times of the stages. Specifically, 

i-1 

m i = m 0 " Yj f 1 + k |)^l + m H f] i = 1, . . . , N (26a) 

£=1 

i-1 

m i = m 0‘ 2 t (1 + k A T l + m H A ~ Vi i = 1, . . • , N (26b) 

i=l 

where hiq is the lift-off mass and m^ ^ and are the fixed hardware mass and 
propellant tank fraction, respectively, for each stage. The summations in equations (26) 
(and all other similar summations) are defined to be zero whenever the lower summation 
limit exceeds the upper limit. The superscripts 0 and f refer to conditions at the 
beginning and the end of each stage, respectively. 

In order to calculate the final payload, the velocity impulse discussed earlier must 
first be considered. This velocity impulse is added after the desired orbit is achieved, 

f H +1-1 

using the N in stage for propulsion. After this maneuver is completed, the N tn stage 
hardware (fixed and variable) is jettisoned. The propellant required to achieve the 
velocity impulse is 
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. m t (, - AV I V t n) 

Anij = m N l 1 - e > J 

and the final payload is 

m PL = m N “ ^ + k N^ Am I “ k N^N T N " m H, N 
The payload can be expressed as a function of burning times, by use of equation (26b). 


m 


PL = y 



+ \ + m 


H,£ 


] - ^N t N f ‘ k N% T N " m H,N ^ 


where y has been defined as 


= -k N + (1 + k N )e 


-AV n 


V t n 


(28) 


0 f 

The variations dm”, dm:, and dmp L can now be expressed in terms of the variations 
in burning times by differentiating equations (26a), (26b), and (27). 


i-1 


dm? = - Yj (! + k £ )% dj i 
£=1 

i-1 

= - Z (1 + k A dT £ -^i dT i 
£=1 


i = 1, . . . , N 


(29a) 


(29b) 


N-l 

dm PL = + k^)^ d-r^ - (y + kjq-)^^ (29 c) 

4=1 


By use of equations (29a) and (29b), the terms in equation (19) involving dn^ can be 
simplified and expressed in terms of the variations dr^: 
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N 



i=2 



- E 4H - < 0 

i=2 


, ,N j f ,1.0 
+ x 5 dm N - X 5 dm 2 


N-l 


- E 4vi dT i 

i=2 


N-l 


4 E (1 + dT i - 4% dT N + 4 1 + k i»i dT i 

i=l 


(30) 


where 


^5 - (^5) 

D 0 t=t. 

The continuity of Ag (from eqs. (16)) is used in equation (30). Also, the time variation 
terms in equation (19) can be expressed in terms of variations in the stage burning times: 

N t N N 

E < c i <*> 1 - E c i< dt i - dt i-i> = E c i dT i (31) 

i=2 l i-l i-2 i=2 

Using the definition of g (eq. (2)) and equation (29c) yields the variation of g: 


N-l 

dg = y Yj ^ + k i^i dT i + ( y + k N^N dr N ( 32 ) 

i=l 

Since a suitable form for all the terms in equation (19) has been obtained, equations (25), 
(30), (31), and (32) are now substituted into equation (19), and the following result is 
obtained : 
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N 


dJ 


-ZwZf 

i=2 i=l x 


X, — + rX 0 — + X„ iL + \ A lSL\ d v . 


1 i A o ' r A S — + a 4 — 1 - | -•/, 

dq- dq. 3 dq. 4 d q-J 3 

1=1 3 3 3 j/ t=t N 


(x 1 ^ + rX ? ^ + X„^- + X 4 ^ da 
\ 1 0a * da 6 da % 9a/ t=t 


L *L + rx ? i" + X, i£-+ x 4 -i£\ 

\ 9t 1 9t 1 9 ^1 dT l) t=t 


dT, 


N-l N-l 

+ E x 5 k i' 3 i dT i - X 5 E <' + k i>' ! i dT i - x 5% dT : 


i=2 


N 


i=l 


N-l 


Xg(l + kj)^ dr i + y (! + dT i + (y + k N )% dT N = 0 (33) 

i=l 


Combining terms yields 
N-l 


dJ 


‘Yj < 1 + k i»i( x - X 5) +C i + X 5 k i' 5 i 


i=2 


dT i 


+ 

(l + k 1 )^ 1 (y-xj +^)-x 1 ^L 

- rk 2 ! w 

-x 3 9r 

- x 4 


L ' / 9tj 

9t 1 

3 9ti 

"lj 


dT, 


N 


(y + k N )/3 N + c N - /3 N x 5 


r N - (* 


dT M -(x,3lL + rx 2 i" + x,2£. 

9a 9a 6 da 


♦ x 4 i£\ 


do? 


S(> 


0U 


rx, ^-+ rx 0 + X n — + X. 

9?7- 

j=l ' 3 


1— ■ +rX 2^ +x 3 7T +X 4^ d ^r° 


9r?j 9??j 9 tj. 


E ) 

^ t=t N 


(34) 
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The coefficient of dr. in equation (34) contains C. and x?. Both of these can be elimi- 

1 N 1 D 

nated by using the following identity for Xg : 


N 


N 


4=i+l XI r/ fi=i+l 




‘ + Z ‘ X 5" 1 ) + 2 ( X 5 ' "s' 1 ) + 4 


—‘Vi 

£=i+l x * 


£=i+l 


i-i+1 

%=0 


(35) 


Equation (7e) shows that x & = 0 for coast phases, so that for = 0, 


4 = 4 _1 


(36) 


It is convenient to define 


, C. 

S = — - - Xg 

h 


s° = — - 4’ 1 
h 


s\ = 0 


S? = 0 


Pi* 0 
Pi * 0 

^i = o 
h- 0 


y i = 2, . . . , N 


3u 


,f _ 1 


Si - (X, — + rx„i^-+ X„ — + x 4 

* 0TJ 6 0TJ 9 Tj y 


(37a) 

(37b) 

(37c) 

(37d) 

(37e) 


With these definitions, equation (35) becomes 
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(38) 


X N - 
A 5 - 


N 

- E 

SL= i+1 


( s \ - s °) 


+ Xi 


By use of equations (37) and (38), equation (34) becomes (if it is assumed that the first 
and last stages are powered): 



The variations in equation (39) are all independent, so that the form of equation (20) has 
been achieved, with 


G(r i ) = (1 + k^jSj 


r + 


E 


L=i+1 


+ ¥i 


^ * 0; i = 1, . . N - 1 (40a) 


0(7^ = C i /3j = 0; i = 2, . . . , N - 1 


G(Tjq-) _ (y + + 


•N^N T F N N 


G(a) = - (A, ^L+rA 9 ^ + A„^- + X 4 ^ 

\ 1 9a z 9a ^ 9a * 9a/ t=t 


(40b) 

(40c) 

(40d) 
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G(„,)=(x. *L + r x.iiUx-^+X.M j = l, (40e) 

1 2 an, 3 s ij Sit 

i- l N 


From equations (11) and (37), 


t=t. 

— j3.^0; i=2 N (41a) 


2 

2 

t=t i 1 

1 — /S.*0; i=2, . . N (41b) 


h = 0; *!_! < t < *1 (41c) 


Since equations (41) do not contain C or X 5 , equation (7e) need not be integrated to 
evaluate equations (40), as indicated earlier. 




Boundary Value Problem 

The determination of an optimum trajectory requires the simultaneous integration of 
the equations of motion (eqs. (3a) to (3d)) and the Euler -Lagrange equations (eqs. (7a) to 
(7d)). A set of initial conditions (state variables and Lagrange multipliers) and staging 
times is required in order to start the integration and to specify the trajectory uniquely. 

The trajectory thus generated must satisfy N + 5 final conditions, corresponding to 
the N + 5 independent problem variables in equation (39). For specified variables x^, 
the final conditions have the form 


x. = x- A 

1 J,d 


( 42 ) 


where the subscript d indicates the desired final value . For unspecified variables, 
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equations (40) supply auxiliary final conditions with the form 


G( X j) = 0 


( 43 ) 


Some of equations (42) are easily satisfied; for example, a specified burning time 
for any stage can be achieved simply by terminating that stage at the proper time during 
the integration. 

An iteration is required in order to satisfy the nontrivial final conditions, and 
variable initial conditions (equal to the number of final conditions) must be available. 

For the present problem, the initial state conditions cannot be varied independently, 
since these variables are determined by the choice of booster burning time and kick angle. 
The Lagrange multipliers (X^, i = 1, . . .,4), however, can be varied independently . 

In addition, the burning times of all stages being optimized and the booster kick angle 
(if it is being optimized) are available as variable initial conditions. 

The initial and final conditions for the two -point boundary value problem can be 
listed as follows: 


Unknown initial condition Desired final condition "1 


X 

X 

X 

X 


1 

2 

3 

4 


a 


¥‘n> = ■’j.dl 


or > j = 1, . . . , 4 


G(hj) = 0 J 

G(a ) = 0 


1 




(44) 


Equation (44) contains K + 5 unknown initial conditions (for K optimized stages, K < N) 
and an equal number of desired final conditions. The size of the iteration loop can be 
reduced by using the fact that equations (7) are homogeneous in the X's. This implies 
that the choice of any one X in equation (44) is arbitrary and serves only as a scale 
factor for the others. The value of this multiplier can be chosen to satisfy one of equa- 
tions (40a) or (40c) appearing on the right side of equation (44). The choice of the arbi- 
trary multiplier is made in appendix B. 
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The iteration size can be further reduced when X 4 = 0 is a required final condition, 
which occurs when the travel angle is unspecified, as shown in appendix C. Since X 4 is 
a constant (eq. (9)) and is continuous across staging (eq. (16)), this final condition can 
be satisfied at t = t^; X 4 is thus removed from both sides of equation (44). 

Equation (40d) can also be evaluated at t = t^, and the booster kick angle optimiza- 
tion can be removed from equation (44). This is accomplished by choosing one X at 
t = t^, such that equation (40d) is satisfied. 

The burning times of optimized stages can sometimes be removed from the iteration, 
along with an equal number of equations (40a), (40b), and (40c). The principle involved 
is similar to that used in removing equation (40d): If all variables in the equation 
G(t£ ) = 0 can be calculated at (or previous to) t = t^, stage SL can be terminated (during 
stage l ) whenever the optimizing equation is satisfied. The number of optimized stages 
which can be thus removed depends on the characteristics of the particular problem. 

Consider first two optimized powered stages and t , i < m, with k m # 0. 
Equation (40a) must be set equal to zero for each stage: 


(1 + k £ ) 



N 


£ 



(1 + O 

v m' 


N 

y+ £ 
i=m+l 





These two equations can be combined to give 


f m-1 

(s* . s?) ♦ Jm- 

1+k ‘ ife/ ; 1+k - 




(45a) 


(45b) 


(46) 


As indicated earlier, equation (45a) can be satisfied by the proper choice of the initial 
scale factor for the X’s. With equation (45a) satisfied, equations (45b) and (46) are 
equivalent; however, equation (46) contains terms which can all be calculated at, or 
prior to, stage m cutoff, and this equation can be used to terminate stage m. Specifi- 
cally, stage m is terminated when 


1 + k 


S = S 1 = 
m m 
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m-1 
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1 + k„ 
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(47) 
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If k m = 0, the term in disappears from equation (46), and this equation can then be 
used to terminate stage m - 1: 



m-1 

z 

i=!+l 


m-2 



i=i+l 


1 + k„ 


m > i + 1 




l + k £ 


m = i + 1 


(48a) 


(48b) 


The terms on the left side of equations (48) are evaluated at stage m-1 cutoff, and the 
terms on the right side of (48a) are evaluated during previous stages. 

For the special case m = N, * 0, equation (40c) 

(y + kjyj)/3j^ + = 0 (49) 


is used instead of equation (45b). For this case, the scale factor is calculated from 

-P 

equation (45a), after is first eliminated by use of equation (49). This scale factor 
is then applied to equation (49) to give 



y + k 
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} + kg. 
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Equation (50) is used to terminate stage N. 

If m = 2 and SL = 1, equation (48b) becomes 


S 


0 

2 ‘ 



1 + kj 


= 0 


(50) 


(51) 


Equation (51) could be used to terminate the booster stage optimally. An alternative 
procedure, however, is to terminate the booster at a prespecified time and to choose 
the initial value of one of the Lagrange multipliers such that equation (51) is satisfied. 
This procedure is similar to that used in eliminating equation (40d) and is computationally 
more convenient, as will be shown in appendix B. 

For coast phases to be optimized, equations (40b) supply criteria for termination of 
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the previous stage, similar to equations (48) and (51). For such cases, stage l - 1 is 
terminated when = 0, where stage 1 is the coasting stage to be optimized. 

In general, use of the preceding equations can supply either zero, one, or two equa- 
tions for termination of a stage. If no equation is supplied, the burning time of that stage 
must be used as an initial condition in the external iteration. If one or two equations for 
termination are available, the stage is terminated by use of one of the equations, and 
the additional equation (if available) is used as a desired final condition in the external 
iteration. 

One other simplification of the two -point boundary value problem is possible when 
stage N is being optimized but cannot be terminated by use of equation (50) (this occurs 
when kjj = 0). For this case, stage N can be terminated when one of the desired orbit 
conditions (appendix C) is satisfied; and one of the orbit conditions are thus elimi- 
nated from equation (44). It is important, however, to choose an orbit condition which 
increases or decreases monotonically during flight, so that the desired value will always 
be achieved. Energy and velocity are two examples of such terminating variables. 


Choice of Initial Conditions 

In order to facilitate convergence of the two-point boundary value problem, it is 
desirable to start the iteration process with initial guesses which are as close as possible 
to the converged values. It is difficult to guess values for the Lagrange multipliers. It 
is possible, however, to use the thrust direction ip and its derivative ip as variable 
initial conditions, and to calculate the initial values of two of the multipliers from these 
variables. The conversion equations are derived in appendix B: 

Xj = Xq sin ip (B4) 

X 2 = Xq cos ip (B5) 

Xq = — - — few - ip - — sin ip cos ip\ + — tan ip ip ^ — (B6) 

15 cos ip \ r / r 2 

where Xq is the arbitrary (positive) scale factor. If either the booster burning time or 
kick angle is being optimized internally, the appropriate optimizing equation is used to 
calculate ip rather than one of the multipliers. When both the booster burning time and 
kick angle are optimized internally, ip and are calculated from the optimizing equa- 
tions. An iteration is required to calculate these variables, and the equations and pro- 
cedures are presented in appendix B. 
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Examples 


The equations and procedures derived in the preceding sections can best be illustrated 
with the aid of examples: 

(1) A three-stage vehicle, flown to circular orbit at a specified radius with the travel 
angle unspecified, where all three stages and the booster kick angle are available for 
optimization and k = 0 for stages 2 and 3. 

Equations (40a) and (40c) must be satisfied for the optimization of the propellant 
loadings: 


3 1 

(1 + kj)^ y + 2 ( S i - S?) + ^s[ = o 

L i=2 X ' J 

(52a) 

(l + k 2 W 2 (y + S £ 3 .S») + ^S [ 2 = 0 

(52b) 

(r + kg)/3g + /3gSg = 0 

(52c) 


Combining these equations as in equation (48b) gives (if kg = kg = 0): 

S f 

n ®i 

S« — = 0 (53a) 

^ 1 + k 1 

S® - S* = 0 (53b) 

Equations (53) are used in place of (52b) and (52c). Since the booster kick angle is also 
being optimized, equations (40d) and (53a) are used to calculate $ and r^, as shown in 
appendix B. Equation (53b) is used to terminate stage 2 internally, and equation (52a) is 
satisfied by the proper choice of the scale factor \q. Since no optimizing equation is 
available for terminating the final stage, this stage is terminated when the desired final 
velocity is achieved. 

The final orbit requirements for this problem are given in appendix C by equa- 
tion (C3). Since can be removed from the iteration and the desired velocity has 
been used to terminate stage 3, only two initial and final conditions remain for the exter- 
nal iteration: 
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Initial 

conditions 


Final 

conditions 


I (54) 

0 r = r d f 

a u = 0 

(2) Same as example (1) but with a fixed second-stage burning time (propellant 
loading): For this case, equations (52a) and (52c) must be satisfied, and these equations 
are combined as in equation (48a) to give 



(55) 


Equation (55) cannot be used to terminate stage 2, since this stage must be terminated 
when t - t^ = T 2 » Also, equation (53a) is not available to terminate stage 1, so that 
must be placed in the iteration. Equation (40d) is used to calculate ip. As in the pre- 
vious example, equation (52a) is satisfied by the choice of scale factor, and stage 3 is 
terminated on final velocity. The iteration variables for this example are 


Initial conditions 


a 


$ 


Final conditions 
u = 0 

r = r , 


of q° _ o° 
b 2 " b 3 ~ b 2 " 




1 + k 


1 J 


(56) 


If kg + 0 in this example, equation (50) would be used instead of equation (55): 

f y + k„/s^ q £ q 

4 = --7t^ 7-V- s 3 + s 2- s 2 


[1 + k 


1 


(57) 


Equation (57) would be used (instead of the final velocity) to terminate stage 3, and the 
final conditions become 
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u = 0 


r " r d (58) 


cu = 


while the initial conditions remain the same as in equation (56). 



PROCEDURE AND IMPLEMENTATION 

In order to obtain numerical results, the equations derived in the preceding sections 
were programed for solution on an IBM 7094 computer. Equations (3) and (7) are inte- 
grated by using a variable step-size Runge-Kutta integration scheme with an error con- 
trol to minimize truncation effects (ref. 18). 

The two-point boundary value problem of equation (44) is overcome by using a multi- 
variable Newton -Raph son iteration scheme. With this method, changes in final condi- 
tions are assumed to be related to changes in initial conditions by first-order, finite- 
difference equations: 


6Y = M 6X 


(59) 


where SX and 6Y are n-vectors (with an n x n iteration assumed) denoting differences 
in initial and final conditions, and M is an n x n matrix of partial derivatives of final 
conditions with respect to initial conditions: 


3Y i 

M * = ^; 


A Yj 


- k AX k 


(60) 


The matrix M is obtained by integrating a reference trajectory and n independent 
perturbed trajectories, so that the equation 

AY = M AX (61) 

is obtained where AX is an n x n matrix of differences in initial conditions such that 
AX-, is the difference of the j n initial condition (from its value on the reference tra- 
jectory) on the k perturbed trajectory, and AY is an equivalent matrix of differences 
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in final conditions. Since the first guesses of initial conditions will not, in general, lead 
to a converged solution, the guesses are improved by using the equation 

6Y r = M SX r (62) 

where the subscript r indicates differences between reference and desired values. The 
predicted changes in initial conditions 6X p are obtained by combining equations (61) 
and (62): 


6X r = AX AY" 1 6Y r (63) 

Equation (63) is used to predict the initial conditions for each iteration cycle, and the 
iteration proceeds until the desired final conditions are satisfied within some prespecified 
tolerance. 


Criteria for Comparison and Termination 

Since the final payload is to be optimized, it is desirable that the payload error, as 
well as the errors in final state conditions, be within tolerance before the iteration is 
terminated. A measure of the error in payload is obtained from the following equation: 



Equation (64) is also used to determine whether or not the iteration is converging. If the 
value of 6m on the new reference trajectory is larger than the 6m on the previous 
reference trajectory, it is likely that the magnitude of 6X r is so large that the higher 
order terms ignored in equation (59) have become important. In such cases, a new tra- 
jectory is flown with a reduced value of 6X r , and this procedure is repeated until the 
value of 6m is less than that of the previous reference trajectory. When 6m has de- 
creased, the associated trajectory is used as the new reference trajectory, and the iter- 
ation proceeds. 


Booster Table 

In many trajectory studies, the booster stage remains fixed while other parameters 
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(such as upper -stage thrust levels and orbit parameters) are allowed to vary. Because of 
the large number of cases that can be run with the same booster stage, it is convenient to 
integrate a family of boost trajectories, and to store the burnout points in table form. 

The required booster burnout conditions can then be obtained from the table (by use of 
some interpolation scheme), rather than from an integration of new booster trajectories 
for each case. With lift-off thrust-to -weight ratio fixed, for example, two booster de- 
grees of freedom are available: kick angle and burning time. The booster table is thus 
a two parameter family of burnout points (radius, velocity, flight path angle, and weight). 

For a problem involving a new booster stage, three booster trajectories are initially 
flown at different kick angles in the estimated region of interest. For each of these tra- 
jectories, the state variables are stored at fixed time intervals, ranging from minimum 
to maximum propellant loading. Other values of kick angle are flown and added to the 
table as required. During each computer run, the new booster data developed is stored 
and added to the table, so that after a series of computer runs using the same booster, a 
comprehensive table of booster burnout points is available. 

A two-dimensional second-order interpolation scheme is used to determine the burn- 
out conditions for intermediate values of kick angle and burning time. The partial de- 
rivatives required in equations (40a) and (40d) are also determined, by differentiating the 
second-order interpolation polynomial. The interpolation accuracy depends on the table 
spacing as well as on the amount and form of variation between table points. In order to 
cover the region of interest with a minimum of trajectories, the initial kick angle spacing 
is large. On the other hand, time spacing is small, since only one trajectory is required 
to obtain a complete set of time points for each kick angle. 

In order to determine whether or not the interpolation accuracy is acceptable, an 
additional booster trajectory is flown after convergence is achieved by using the converged 
values of burning time and kick angle. The upper stages are then reflown with the exact 
booster burnout conditions. The payload and burnout conditions from the exact flight are 
then compared with the results of the interpolated converged flight. If agreement is ac- 
ceptable, the interpolated results are retained and the problem is finished. If agreement 
is not acceptable, new booster trajectories are flown and added to the table (in the region 
of convergence) by using valves of kick angle halfway between the previous table points. 
The flight is then reconverged with the modified table. The interpolation accuracy is 
tested as before, and the entire process is repeated (if necessary) until acceptable agree- 
ment is reached. 


Convergence Properties 

When the problem to be solved is not well known to the analyst, it is difficult to obtain 
initial conditions which are close to their converged values. On the other hand, the 
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equations for optimal staging are sometimes very strongly coupled to the initial condi- 
tions. For this reason, it is desirable to optimize the steering program and booster kick 
angle first, before attempting to optimize burning times. 

For a three-stage vehicle to be completely optimized (booster kick angle and all three 
burning times optimized), the procedure is as follows: First, the booster and second- 
stage burning times are held fixed at initial estimated values while the booster kick angle 
and third-stage burning time are optimized. The converged values of kick angle, pitch 
rate, and third-stage burning time are then used as first guesses, and the second- and 
third-stage burning times are optimized in the second case. Finally, all three stages 
(and the booster kick angle) are optimized in the third case. 

This procedure has been used with good results in solving problems which were pre- 
viously unfamiliar (as in the results presented). Of course, when the problem to be 
solved represents a small perturbation from a previously solved problem, complete opti- 
mization can be accomplished in one pass. 


Perturbation Size 

As discussed earlier, a finite difference procedure is used to obtain the partial de- 
rivatives of final conditions with respect to initial conditions. The perturbation sizes 
used for the initial conditions must be carefully selected: if the perturbations are too 
small, round-off and truncation errors inherent in the integration scheme cause large 
errors in the final differences; if the perturbations are too large, nonlinear effects be- 
come important, and incorrect local derivatives are obtained. The difficulty in choosing 
the proper perturbation sizes is magnified by the strong coupling of the partial derivatives 
to the optimal staging equations. 

In order to obtain accurate local partial derivatives, an automatic perturbation size 
checking scheme was implemented in the program. With this scheme, the perturbation 
sizes are adjusted (when required) so that the payload error difference (eq. (64)) between 
reference and perturbed flights is between fixed limits. Generally, the perturbation 
sizes must be adjusted whenever the optimization criteria are changed. 


RESULTS 

In order to demonstrate the validity of the equations and the feasibility of the vari- 
ational technique, parametric results are presented and compared with the overall opti- 
mum solutions obtained by using the variational technique. Two- and three-stage launch 
vehicles were optimized by use of the equations developed. The results presented include 
a two-stage vehicle flown to circular orbit and a three-stage vehicle flown to Earth escape 
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through a circular parking orbit. In these results, both fixed and variable hardware 
weights were used, and the propellant loadings and booster kick angle optimized. Para- 
metric results are also presented, which show the variation of payload with propellant 
loadings and kick angle. 


Vehicle Definition 


The vehicle chosen for this study is a hypothetical three-stage launch vehicle con- 
sisting of two chemical stages and one nuclear stage. The assumptions on propulsion and 
weights are listed in table I. The first stage is based on the Saturn SI-C stage, consist- 
ing of five F-l engines using RP-LOX propellants. Each engine delivers 1. 5 million 
pounds of thrust at sea level, with a specific impulse of 264 seconds. These data are 
used for illustrative purposes only and are not necessarily consistent with present SI-C 
values. The second stage consists of one M-l engine using liquid hydrogen and liquid 
oxygen as propellants. This stage operates at a vacuum thrust of 1. 5 million pounds and 
a vacuum specific impulse of 428 seconds. The third stage is a nuclear stage, with esti- 
mated thrust of 250 000 pounds and a specific impulse of 850 seconds. 

As explained earlier, the booster burnout conditions are determined from a table as 
functions of burning time and kick angle. The launch thrust -to -weight ratio is fixed at 
1.25, so that the launch weight is 6 million pounds. The booster stage is flown vertically 
for 15 seconds, at which time the relative velocity vector is instantaneously tipped to the 


TABLE I. - LAUNCH VEHICLE PROPULSION 


AND WEIGHT DATA 



Stage 



First 

Second 

Third 

Thrust, lb 

7. 5xl0 6 (Sea level) 

1. 5xl 0 6 

2. 5x10® 

Specific 

264 (Sea level) 

428 

850 

impulse, sec 

305 (Vacuum) 



Fixed 

245, 000 

70, 000 

35, 000 

hardware 




weight, lb 




Propellant 

0. 030 

0.033 

0.120 

sensitive 




fraction 




Drag reference 

855 



area, sq ft 






Figure 2. - Drag coefficient as function of Mach number. 
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Figure 3. - Payload capability as function of second-stage propellant 
loading. Three stages to Earth escape via 121-nautical-mile park- 
ing orbit; variable tanks; first-stage propellant loading at optimum 
point, 4 028 000 pounds; third-stage propellant loading at opti- 
mum point, 317 000 pounds. 



First-stage propellant loading, lb 

Figure 4. - Payload capability as function of first-stage 
propellant loading. Two stages to 121-nautical-mile 
orbit. 


desired kick angle (measured from 
the horizontal) and azimuth heading. 
A launch azimuth of 90° is assumed. 
Following the initial vertical rise 
and pitchover, the booster is flown 
with a zero angle of attack steering 
program. 

The booster trajectory simula- 
tion includes a detailed oblate Earth 
gravity model (ref. 19). The Air 
Research and Development Command 
atmosphere model (ref. 20) is as- 
sumed. Booster drag is calculated 
as a function of Mach number by 
using the drag coefficient curve pre- 
sented in figure 2 (p. 29). 

The second and third stages use 
the calculus of variations steering 
program, as explained earlier. A 
spherical Earth with no atmospheric 
forces is assumed for these stages. 


Numerical Results 

A typical propellant tank sizing 
study was conducted with the vehicle 
defined in table I. In this study, the 
stages are assumed to have variable 
tank size, and the optimum propellant 
capacities are determined by flying 
possible missions of interest. 

The first mission is flown to 
Earth escape energy with three 
stages. This mission is typical of 
lunar and planetary probes and or- 
biters. A parking orbit ascent mode 
is used, wherein all three stages are 
used to enter a parking orbit of a 
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I 



89.30 89.35 89.40 89.45 89.50 89.55 89.60 


Booster kick angle, deg 

Figure 5. - Payload capability as function of booster kick angle. Pro- 
pellant loadings optimized; two stages to 121-nautical-mile orbit; 
first-stage propellant loading at optimum point, 3 980 000 pounds; 
second-stage propellant loading at optimum point, 1 276 000 pounds. 

121 -nautical-mile altitude (similar to the Apollo mission), after which the third stage is 
burned to escape. The second burning of the third stage is assumed to be impulsive, 
with a velocity increment of 10 563 feet per second. The maximum payload capability for 
this mission is obtained by optimizing the booster kick angle as well as the propellant 
loadings of the three stages. By use of the equations derived earlier, maximum payload 
can be obtained with a single (converged) solution, represented by the optimum point in 
figure 3. Actually, since the hypothetical vehicle used for illustration is not typical of 
multistage chemical vehicles, good estimates of initial conditions (in this case, ot , t^, 
and ip) were not known. The three-step convergence procedure discussed earlier was 
thus used. 

With the parametric procedure, the maximum payload is obtained as the envelope of 
payload points with all possible combinations of propellant loadings and booster kick 
angle. The increased effort (number of solutions) required for this procedure is obvious. 
In addition to the parametric curves shown in figure 3, each point on the curves had to be 
obtained at an optimum kick angle, which required an additional family of curves (not 
shown in the figure), with kick angle as the independent parameter. As can be seen from 
figure 3, the payload capability and optimized propellant loading obtained from the varia- 
tional procedure agree with the parametric results. 

Another mission of interest is two stages flown to a circular orbit, with the use of 
the first and second stages from table I (p. 29). Results for this mission are presented 
in figures 4 and 5. In figure 4, with the use of the parametric procedure, payload capa- 
bility is presented as a function of first-stage propellant loading for various booster kick- 
angles. The maximum payload for each kick angle is determined from the figure and 
presented as a function of kick angle in figure 5. The case flown by use of the variational 
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Figure 6. - Payload capability as function of second-stage propellant 
loading. Three stages to Earth escape via 121- nautical-mile park- 
ing orbit; fixed tanks; first-stage propellant loading at optimum 
point, 3 964 000 pounds; third-stage propellant loading at optimum 
point, 310 000 pounds. 


technique is represented by the opti- 
mum point in figure 5, and the pay- 
load obtained is in agreement with 
the envelope, as before. 

Frequently, the propellant tanks 
of the various stages are sized for 
one mission and fixed at these values 
for all other missions. Propellant 
loadings can still be optimized in 
such cases but with the restriction 
that the propellant loadings cannot 
exceed the propellant capacities of 
the tanks. 


In figure 6, the three-stage Earth escape mission is reoptimized with fixed tanks for 
the first and second stages. The tank weights and propellant capacities used are based 
on optimum values for the two-stage orbit mission (figs. 4 and 5). Since the maximum 
propellant capacities of the first and second stages were not exceeded in the optimum 
case, the propellant loadings of these stages were off-loaded to the optimum valves. The 
parametric results and the variational point are compared (as in fig. 3, p. 30), and the 
results are in agreement. 


CONCLUDING REMARKS 

A technique is presented which allows simultaneous optimization of the thrust direc- 
tion profile and vehicle control parameters for multistage launch vehicles. The agree- 
ment of parametric and optimum results presented in figures 3 to 6 (pp. 30 to 32) demon- 
strates the correctness of the optimizing equations. 

The amount of effort (and computer time) saved by the variational technique in de- 
termining the maximum payload capability can be seen by referring to the figures. In 
figure 3, for example, approximately 80 parametric data points are required to optimize 
the three-stage propellant loadings and booster kick angle (not shown in the figure). When 
good initial guesses are available, the computer time required to obtain the overall opti- 
mum solution is about the same as for each parametric solution. The time saving is 
therefore proportional to the number of parametric cases required for complete optimi- 
zation. Even when good initial conditions are not available (this was the case in fig. 3), 
three cases are usually sufficient to obtain the optimum solution (as described earlier), 
and the time saving is still substantial. 

The technique presented is, of course, not limited to the parameters (propellant 
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loadings and booster kick angle) used herein. Other possible parameters (for launch ve- 
hicle optimization problems) are launch thrust -to -weight ratio and thrust levels for the 
various stages. In addition, the technique should be applicable to other types of optimi- 
zation problems. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, October 8, 1965. 
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APPENDIX A 


SYMBOLS 


a functions defined in eqs. (BIO), 
ft/ sec 2 

b functions defined in eqs. (B13), 
sq ft/(sec 2 )(rad) 

C constant of integration 

E Weierstrass excess function 

e eccentricity 
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j=l 

f constraint equation 
g aj/ax 

g function of initial and final conditions 
to be minimized 

2 

H energy per unit mass, sq ft/ sec 

J functional to be minimized by varia- 
tional methods 

k propellant sensitive mass fraction 

M matrix relating changes in final and 
initial conditions 

m mass, slugs 

p semilatus rectum, ft 

r radius, ft 

S functions defined in eqs. (37) 

T thrust, lb 

t time, sec 

u radial velocity, ft/ sec 
v velocity, ft/ sec 

x problem variable 


a booster kick angle, rad 
/3 mass flow rate, slug/sec 
T flight path angle, rad 
y function defined in eq. (28) 
e argument of perigee, rad 
77 generalized state variable 
6 true anomaly, rad 
X Lagrange multiplier 
Aq arbitrary scale factor 

2 

M Earth force constant, cuft/sec 
t burning time, sec 
c p polar angle, rad 

ip thrust direction, rad 

o> angular velocity, rad/ sec 

Subscripts: 
d desired value 

H fixed hardware 

I impulsive 

i stage number 

j variable number 

k variable number 

i stage number 

N last stage 

PL payload 

p propellant 

r difference between reference and 
desired values 
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s 

structure 

derivative with respect to time 

0 

initial 

* finite, but allowable variation of the 

Superscripts: 

variable 

f 

end of stage 

+ after staging 

i 

refers to t = t. 

before staging 

0 

beginning of stage 

— vector 
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APPENDIX B 


CONVERSION EQUATIONS 
Calculation of Lagrange Multipliers 

In order to avoid the difficulty associated with guessing initial values of the Lagrange 
multipliers, equations are derived which express two of the multipliers in terms of the 
pitch attitude ip and rate ip. 

Equations (8a), (8bl) and (8cl) give ip in terms of the Lagrange multipliers: 

X 1 

tan ip = — (8a) 

X 2 

X I 

sin ip - - (8b 1) 

V x I + 4 

Xn 

cos i[) = — , (8cl) 

V*i + 4 


Differentiating equation (8a) yields the pitch rate 


'/ . X 2 X 1 ' X 1 X 2 
ip sec ip = 


which becomes, with the use of equations (8cl), (7a), and (7b), 

(v 3 * ; x i x 2 - ^ 


if . i = 2(0 - 


4 +X 2 


(Bl) 


(B2) 


The discussion following equation (44) shows that one of the multipliers can be picked 
arbitrarily. It is convenient to pick Xq as the scale factor, where 
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X 0 ~ V x 1 + x 2 ^ 0 


(B3) 

With this choice, 







Xf = Xq sin Ip 


(B4) 



X£ = Xq cos ip 


(B5) 

Substituting equations 

(B3), 

(B4), and (B5) into (B2) and solving for 

Xg gives 


^3 = 

X 0 

(2 co - ip sin ip cos ; p\ + — tan ip 


(B6) 


COS Ip 

\ r Jr 

2 


The multipliers X^ 9 ) 

\_ 2 , and Xg can now be calculated from ip, ip, 

X 4 , and Xq 



by using equations (B4), (B5), and (B6). 


Internal Optimization 

The analysis of the boundary value problem showed that the booster kick angle and, 
in certain cases, the booster burning time can be optimized without iteration. In such 
cases, the optimizing equation can be used to calculate ip or ip; one set of initial and 
final conditions is thus eliminated from the iteration. 

The equation for optimizing the booster kick angle is 


3u . 3co , . 3r 3 cp n 

X-, + rX 9 — + Xq + Xa — = U 

1 3a z 3a 6 da * da 


(40d) 


or, in terms of ip and ip , 


3u 


3 co 


X 0 sin ip — + rx 0 cos ip — + 


30! 


3a 


[cos ip 


0 / • u \ ^4 

^2 co - ip sin ip cos ip) + — 


tan ip 


1L +X 4 ^ = o 

da da 


(B7) 
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Equation (B7) can be solved for ip: 


ip = 2co - — sin ip cos ip + — s * n ^ + cos ^ /sin ip — + r cos ip — + -i (B8) 

r X Q r 3r_ l da 3 a Xq 3 a) 

da 

Use of equation (B8) to calculate ip, however, leads to serious convergence problems, 
since this equation is an extremely nonlinear function of ip. Instead, ip is used as a 
variable initial condition in the external iteration and equation (40d) is used to calculate 
ip. Although a closed form solution for is not available (unless X^ = 0), an iteration 
can be performed which converges rapidly. The procedure is as follows: 

(1) Guess ip. 

(2) Calculate x^, Xg, and Xg from equations (B4), (B5), and (B6). 

(3) Calculate A = X 1 — + rx 9 — + X„ — + X 4 ^ . 

1 So; * da 0 da * da 

(4) Iterate steps (1) to (3) to obtain A = 0. 

The booster burning time can be optimized internally when 

(1) Stage 2 is a powered stage (j3g, Tg # 0) also being optimized and kg = 0 

(2) Stage 2 is an optimized coast phase 

The optimizing equations for these two cases are 

S? - = 0 (51) 

1 1 + k x 


and 


C 2 = 0 (40b) 

If equations (37e), (41b), and (41c) are used and it is noted that 3r/3r^ = u and 
dcp/dr i = co, equations (51) and (40b) become 

a^X^ + agXg + a 3 Xg + a^X 4 = u^Xq (B9) 


where 
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(BIO) 


Equation (B9) can be expressed in terms of ip and ip as in equation (B8): 


ip = 2 u sin ip cos ip + cos t /a x sin ip + a 2 cos ip + — a 4 - a 5 | (Bll) 


X 0 r 


As in equation (B8), however, ip is a very nonlinear function of ip, and equation (B9) 
must instead be solved for ip in terms of ip. Again, a closed form solution is not avail- 
able, and the iteration proceeds as previously described, except that 


A - ajXj + a 2 X 2 + a 3 X 3 + a 4 X 4 - a 5 X Q 


for this case. 

If the booster burning time and kick angle are both being optimized internally, equa- 
tions (40d) and (B9) are available for the calculation of any two initial conditions. It is 
convenient to use these equations to calculate ip and Tj and to leave ip and a as vari- 
able initial conditions in the external iteration. 

Equations (B9) and (40d) are combined to eliminate Xg, after which the resulting 
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equation is used with equation (B3) to solve for X-^ and X 2 : 



(B12a) 


(B12b) 


where 

b 1 = a< — - a, (B13a) 

1 1 3a 3 da 

b 0 = a, — - a„r — (B13b) 

6 * da 6 da 

N= fa^-a 4 — V 4 +a 5 X 0 ^- (B13c) 

3 \ 3 Bar 4 * da) 4 3 U 0O! 

The choice of sign in front of the radical in equations (B12) was made to obtain a thrust 
direction between 0° and 90°. Substituting equations (B4) and (B5) into (B6) gives 

Xn = — (2u> - ip) - — X-, + — (B14) 

3 X 2 r rx 2 

An iteration is required to calculate r^, which proceeds as follows: 

(1) Guess Ty 

(2) Calculate Aj, X 2 , and Ag from equations (B12) and (B14). 

(3) Calculate A = A, — + rX 9 — + X, — + X 4 

1 da 3 da 3 da da 

(4) Iterate steps (1) to (3) to obtain A = 0. 

Equations (B12) and (40d) are equivalent to equations (B3), (B9), and (40d), and thus, the 
optimum values of ip and are obtained when the iteration converges. The value of 

ip is calculated from equations (B12), by use of (B4) and (B5): 
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(B15a) 


sin ip = 


b l b 3 



V( b 1 + b 2>0 

o( b l + b 2) 



COS Ip = 


b 2 b 3 +b iVMi!lM 

*o( b I + b !) 



(B15b) 


The value of ip obtained from equations (B15) is independent of the value of Aq chosen, 
provided that A^ is scaled exactly as Aq. 
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APPENDIX C 


ORBIT REQUIREMENTS 

A two-dimensional orbit is completely specified by any four independent orbit param- 
eters. In many cases, however, one or more of these parameters is left unspecified and 
can be used to maximize the payload. For such cases, equations (40e) supply auxiliary 
boundary conditions to be satisfied, in order to guarantee optimum values of the unspeci- 
fied variables . 

A typical example is the requirement of a circular orbit at a specified radius, with 
the orbit injection point unspecified. For this case, the required final conditions can be 
stated as 



(Cl) 


where r, u, and o> (and all other variables in this appendix) should be evaluated at 
t = t^. Equations (40e) are written (with r, u, co, and cp as independent variables) 


G(u) = 

G (u>) = rx 2 
G(r) = X 3 




G(<p) = 


(C2) 


Since <p is not specified, X^ = 0 is the auxiliary boundary equation, and the required 
final conditions are 
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r 


= r ■ 


u = 0 


Cl) = 



> 


(C3) 


X 4 


= 0 


If, in addition to the travel angle, the radius is also unspecified, equation (C2) gives 
Xg = 0 as an additional auxiliary boundary equation, and 

u = 0 ^ 



(C4) 


are the required final conditions. 

Frequently, the orbit is specified in terms of energy and flight path angle. For such 
cases, H, r, r, and cp can be used as independent variables, where 


H = — (u 2 + r 2 a> 2 ) - — 


and 


(C5) 


tan r = — 
cur 
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Solving for u and oj yields 


u = y 2 (H + iij sin r 


co =1^2^ + -^cos r 




(C6) 


If the indicated partial derivatives are taken, equations (40e) become 


g(h)= & xi+ :h 


G(r) = (corX^ - uAg) 


G(r) 


im 

2 2 
v r 


(- 5 ) 




x 2 + x 3 


G ((p) = *-4 


(C7) 


If energy and flight path angle are specified (radius and travel angle optimized) the re- 
quired final conditions are 


- (u 2 + w 2 r 2 ) - — = H, 
2 r a 


tan-ViLU 

\ 0?r 


2 2 
v r 


:)- r d 

(•*3- 


>- 


(C8) 


MU Xj - (u + -^]X 2 + x 3 “ 0 


X 4 = 0 


J 


or, if energy is specified and all other variables are optimized, the required final con- 
ditions become 
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1 / 2 2 2\ u 

- (u + r a) ) - — = H j 

2 r d 


wrXj - uXg = 0 


jLtU 

2 2 
v r 


+ ~^2 j* 2 + X 3 " ° 


x 4 = 0 


(C9) 


The last three of equations (C9) are equivalent to 

*\ 

ip = r 

i = r 


x 4 = o 


(CIO) 


The case of an elliptic orbit is discussed in reference 8, in which the orbit elements 
e, p, 6 and e are chosen as the independent variables: 




e sin 9 


■V? 


co = /— (1 + e cos qY 


(Cll) 


r = 2 

1 + e cos 8 

cp = 6 + € 


For this case, equations (40e) become 
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2 o 

u , ,2r w cos 0 . r cos 0 . 

G(e) = — A-. + Ag Ag 

e p P 


p J 2p 1 2p J 


0(9) = g - r2eal i 1 ” ? X 2 X f2e Bin 8 X, x X 4 


tan 0 


G(e) = X 4 


(C12) 


Frequently, the eccentricity e and the semilatus rectum p of the ellipse are specified, 
while the true anomaly 0 and argument of perigee e are left open for optimization. For 
this case, the required final conditions are 


4 2 
r c o 

P = = P d 

JU 


6=11 + 


2 2 2 o u 4 2 

u + u) r - 2 st r cu 

r 


= e ■ 




U ^1 2r^eco sin 0 r^e sin 0 . _ n 

A 2 + A 3 - u 


tan 0 


a 4 = 0 


(C13) 


46 



REFERENCES 


1. Goldsmith, M. : On the Optimization of Two -Stage Rockets. Jet Prop. , vol. 27, 

no. 4, Apr. 1957, pp. 415-416. 

2. Schurmann, Ernest E.H. : Optimum Staging Technique for Multistaged Rocket Ve- 

hicles. Jet Prop., vol. 27, no. 8, Aug. 1957, pp. 863-865. 

3. Subotowicz, M. : The Optimization of the N-Step Rocket with Different Construction 

Parameters and Propellant Specific Impulses in Each Stage. Jet Prop. , vol. 28, 
no. 7, July 1958, pp. 460-463. 

4. Hall, H.H.; and Zambelli, E.D. : On the Optimization of Multistage Rockets. Jet 

Prop., vol. 28, no. 7, July 1958, pp. 463-465. 

5. Weisbord, L.: Optimum Staging Techniques. Jet Prop., vol. 29, no. 6, June 1959, 

pp. 445-446. 

6. Cobb, Edgar R. : Optimum Staging Technique to Maximize Payload Total Energy. 

ARS J., vol. 31, no. 3, Mar. 1961, pp. 342-344. 

7. Coleman, John J. : Optimum Stage- Weight Distribution of Multistage Rockets. ARS 

J., vol. 31, no. 2, Feb. 1961, pp. 259-261. 

8. Zimmerman, Arthur V. ; MacKay, John S. ; and Rossa, Leonard G. : Optimum Low- 

Acceleration Trajectories for Interplanetary Transfers. NASA TN D-1456, 1963. 

9. Melbourne, William G. ; and Sauer, Carl G. , Jr. : Optimum Thrust Programs for 

Power -Limited Propulsion Systems. Rept. No. TR 32-118, Jet Prop. Lab., 
C.I.T., June 15, 1961. 

10. MacKay, John S.; and Rossa, Leonard G. : A Variational Method for the Optimiza- 

tion of Interplanetary Round-Trip Trajectories. NASA TN D-1660, 1963. 

11. Jurovics, Stephen: Optimum Steering Program for the Entry of a Multistage Vehicle 

into a Circular Orbit. ARS J. , vol. 31, no. 4, Apr. 1961, pp. 518-522. 

12. Mason, J.D.; Dickerson, W.D.; and Smith, D.B.: A Variational Method for Opti- 

mal Staging. Paper presented at AIAA 2nd Aerospace Sciences Meeting, Paper 
No. 65-62, AIAA, Jan. 1965. 

13. Denbow, C.H. : Generalized Form of the Problem of Bolza. Ph.D. Thesis, Univ. 

of Chicago, 1937. 

14. Hunt, R.W.; and Andrus, J.F.: Optimization of Trajectories having Discontinuous 

State Variables and Intermediate Boundary Conditions. Paper Presented at SIAM 
Meeting, Monterey, (Calif.), Jan. 1964. 


47 



15. Stancil, R.T.; and Kulakowski, L. J. : Rocket Boost Vehicle Mission Optimization. 

ARS J. , vol. 31, no. 7, July 1961, pp. 935-942. 

16. Bliss, G.A. : Lectures on the Calculus of Variations. Univ. Chicago Press, 1946. 

17. Leitmann, G. : On a Class of Variational Problems in Rocket Flight. J. Aero/Space 

Sci., vol. 26, no. 9, Sept. 1959, pp. 586-591. 

18. Strack, William C . ; and Huff, Vearl N. : The N-Body Code - A General Fortran 

Code for the Numerical Solution of Space Mechanics Problems on an IBM 7090 
Computer. NASA TN D-1730, 1963. 

19. Clarke, Victor C. , Jr. : Constants and Related Data for Use in Trajectory Calcula- 

tions as Adopted by the Ad Hoc NASA Standard Constants Committee. NASA 
CR-53921, 1964. 

20. Minzner, R.A.; Champion, K.S.W. ; and Pond, N.L.: The ARDC Model Atmosphere, 

1959. Air Force Surveys in Geophysics No. 115, Rept. No. TR-59-267, Air Force 
Cambridge Res. Center, Aug. 1959. 


48 


NASA- Langley, 1966 E-3154 



"The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof 

—National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and technical information considered 
important, complete, and a lasting contribution to existing knowledge. 

TECHNICAL NOTES: Information less broad in scope but nevertheless 

of importance as a contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: Information receiving limited distri- 

bution because of preliminary data, security classification, or other reasons. 

CONTRACTOR REPORTS: Technical information generated in con- 

nection with a NASA contract or grant and released under NASA auspices. 

TECHNICAL TRANSLATIONS: Information published in a foreign 

language considered to merit NASA distribution in English. 

TECHNICAL REPRINTS: Information derived from NASA activities 

and initially published in the form of journal articles. 

SPECIAL PUBLICATIONS: Information derived from or of value to 

NASA activities but not necessarily reporting the results of individual 
NASA-programmed scientific efforts. Publications include conference 
proceedings, monographs, data compilations, handbooks, sourcebooks, 
and special bibliographies. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION DIVISION 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 



